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Abstract 

Analytical forces have been derived in the La- 
grangian framework for several random phase ap¬ 
proximation (RPA) correlated total energy meth¬ 
ods based on the range separated hybrid (RSH) 
approach, which combines a short-range den¬ 
sity functional approximation for the short-range 
exchange-correlation energy with a Hartree-Fock- 
type long-range exchange and RPA long-range 
correlation. The RPA correlation energy has 
been expressed as a ring coupled cluster doubles 
(rCCD) theory. The resulting analytical gradi¬ 
ents have been implemented and tested for geom¬ 
etry optimization of simple molecules and inter- 
molecular charge transfer complexes, where inter- 
molecular interactions are expected to have a non- 
negligible effect even on geometrical parameters 
of the monomers. 

Introduction 

After having achieved the limit of chemical accu¬ 
racy with the hybrid functionals for a vast majority 
of applications, density functional theory (DFT) 
in its Kohn-Sham, 1 or more precisely general¬ 
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ized Kohn-Sham 2 formulation became the most 
popular electronic structure method of computa¬ 
tional chemistry 3 and computational material sci¬ 
ences. 4 In the mean time, new applications re¬ 
vealed more and more examples where routine 
DFT methods fail and pointed out some inherent 
weaknesses of the usual local and semi-local ap¬ 
proximation of the exchange-correlation function¬ 
als. 5 During the last years, considerable progress 
has been achieved in a better understanding of 
the fundamental reasons of these failures, 6 e.g. 
the lack of a proper description of London dis¬ 
persion forces or the systematic errors due to the 
erroneous treatment of charge and spin localiza¬ 
tion/delocalization. 7 Various methods have been 
proposed to remedy functionals for one or more 
of the above-mentioned problems. Among these 
approaches, methods that use both occupied and 
virtual orbitals, sometimes called methods belong¬ 
ing to the 5th rung of Jacob’s ladder, 8 play a privi¬ 
leged role due their extraordinary flexibility, lim¬ 
ited only by the associated computational costs. 
Simple second order Rayleigh-Schrodinger per¬ 
turbation theory can be an acceptable solution 
for many cases, but one runs into serious dif¬ 
ficulties in small-gap systems. The next level 
of approximation, which has a long-standing tra¬ 
dition in DFT, 9 is the random phase approxi¬ 
mation (RPA). The RPA can be derived either 
in an adiabatic-connection/fluctuation-dissipation- 
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theorem framework 10,11 or as a particular approx¬ 
imation to the coupled cluster doubles (CCD) ex¬ 
pression of the correlation energy in diagrammatic 
perturbation theory. 12 The two types of formula¬ 
tions are intimately related and in some special 
cases they lead to correlation energies that are 
strictly equivalent. 13 The great advantage of RPA 
is that it provides a correlation energy functional, 
which is compatible with Hartree-Fock exchange. 

Random phase approximation (RPA) methods, 
belonging to the 5th rung of Jacob’s ladder 8 of 
DFT approaches, are becoming a practical tool to 
construct correlation energy functionals 10,11,14-32 
not only in computational material sciences but 
also in the quantum chemistry applications. Re¬ 
cent works by Eshuis and Furche, 11,33 Hessel- 
mann, 28 Ren et al. 29 and others 20,21 have demon¬ 
strated that highly efficient RPA implementations 
with favorable scaling properties with the system 
size are conceivable. Thus, it is now becom¬ 
ing possible to study extended systems 28,33 with 
similar computational resources but with a con¬ 
siderably better accuracy than the MP2 method. 
Quite good results could be obtained for reaction 
energies and barriers 34 and especially good per¬ 
formance is expected for systems where van der 
Waals interactions, 35 and in particular London dis¬ 
persion forces, 36 play a crucial role. 

A few shortcomings of the RPA method have 
also been identified since the first numerical stud¬ 
ies on molecular systems 14-16 and solids. 20,21,24,37 
As it has been expected from earlier work of 
Perdew and his colleagues, 38 the RPA describes 
rather poorly the short-range correlation: for this 
reason, e.g. its performance for atomization en¬ 
ergies is inferior to that of a good gradient cor¬ 
rected functional. 14 Another problem is due to 
the slow basis set convergence of the RPA corre¬ 
lation energy, which makes the determination of 
CBS (complete basis set) limit correlation ener¬ 
gies (and energy differences) for larger systems 
prohibitively expensive. 17,39 These problems can 
be fixed to a considerable extent by applying in¬ 
stead of the previously mentioned "full-range" 
RPA techniques, which use Kohn-Sham orbitals 
obtained by standard (usually PBE) exchange- 
correlation functionals, followed by the frozen- 
orbital evaluation of the Hartree-Fock (HF) ex¬ 
change and the RPA correlation, the combina¬ 


tion of a short-range DFA (density functional 
approximation) and a long-range RPA method 
within the range-separated hybrid (RSH) frame¬ 
work. 30-32,40^2 

This family of methodologies consists in sepa¬ 
rating the long- and short-range electron-electron 
interactions in the Hamiltonian. The short range 
interactions are described within a density func¬ 
tional approximation (DFA), using appropriately 
designed short-range LDA or PBE exchange- 
correlation functionals. The long-range exchange 
is taken into account by a nonlocal Hartree- 
Fock operator, while the long-range correlation 
is described in the random phase approxima¬ 
tion. Range-hybrid RPA methods have shown 
quite good performance for intermolecular com¬ 
plexes, especially when certain RPA variants are 
applied for the long-range correlation. 43 In par¬ 
ticular, the stacking interactions, which are due to 
long-range dynamical correlation effects, are well 
described. Considering the nature of the RPA, it 
is expected that its best theoretical performance 
is deployed for long-range dynamical correlations, 
which are physically responsible for London dis¬ 
persion forces. Since the reference determinant 
for the long-range RPA calculations is constructed 
from RSH orbitals, which are optimized using a 
long-range HF exchange potential, spurious delo¬ 
calization effects, common for conventional local 
and semi-local functionals can be avoided from the 
outset. In fact, the presence of the long-range HF- 
exchange reduces considerably the delocalization 
error, which usually deteriorates the description of 
charge transfer (donor-acceptor) complexes. Fur¬ 
thermore, since RSH orbitals fulfill a long-range 
Brillouin-theorem, 44 there is no need for single¬ 
excitation corrections, as proposed by Ren and his 
coworkers. 45 It has been shown that the basis set 
convergence properties of the range-hybrid RPA 
method is much faster than that of the "full-range" 
correlation methods and the basis set superposition 
error has an almost negligible impact on the results 
even at relatively small basis sets. 46 Although ac¬ 
tual computational implementations of RSH+RPA 
approaches are far from being optimal, this class of 
methods remains one of the most promising ways 
to correct a number of notorious shortcomings of 
conventional DFAs. 

We should be aware of the fact that the effi- 
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cient and accurate calculation of total energies is 
not enough to perform realistic modeling work: 
one should definitely go beyond total energies and 
in order to relax atomic positions (geometry op¬ 
timization) or follow the evolution of the sys¬ 
tem by solving the equations of motion for the 
nuclei (Born-Oppenheimer molecular dynamics) 
one needs the corresponding analytical total en¬ 
ergy derivatives (gradients) in each point of the 
potential energy surface. Since the pioneering 
work of Pulay on analytical Hartree-Fock deriva¬ 
tives, 47 which opened the way to a routine ana¬ 
lytical calculation of forces, numerous theoretical 
and computational developments have been pro¬ 
posed, making analytical gradients accessible for 
a vast majority of mainstream electronic structure 
methods. The formal difficulties which appeared 
at the beginning for non-variational energy expres¬ 
sions have been overcome by the Lagrangian for¬ 
mulation of derivative properties, 48,49 permitting 
to have analytical forces in ground and excited 
states. 50 Two very recent publications on RPA an¬ 
alytical gradients, which appeared after complet¬ 
ing the present work, have been based also on 
the Lagrangian method. Analytical gradients of 
post-Hartree-Fock RPA correlation energies have 
been published by Helgaker’s group 51 in a ring 
coupled cluster type formulation of RPA, 12 while 
Furche and his coworkers 52 used a direct RPA cor¬ 
relation energy expression, formulated via the fre¬ 
quency dependent dielectric matrix in the resolu¬ 
tion of identity approximation, using Kohn-Sham 
orbitals. Our approach, developed independently 
from theirs, follows a Lagrangian strategy too. 

Our main objective was to derive RSH+RPA 
gradients, i.e. the analytical first derivative of the 
combination of a short-range DFA and a long- 
range RPA method within the range-separated hy¬ 
brid (RSH) framework. We propose gradient ex¬ 
pressions for different variants of the long-range 
RPA correlation energy, 43,53 without and with ex¬ 
change. Our expressions, at the extreme limits 
(zero and infinity) of the range-separation param¬ 
eter, become identical either with post-Hartree- 
Fock RPA, or with pure DFA analytical gradients. 
The second-order limit of the RSH+RPA correla¬ 
tion energy with exchange is the RSH+MP2. 44 For 
this latter case, more precisely for the long-range 
local MP2 approach 54,55 combined with short- 


range DFT, analytical gradients are already avail¬ 
able, 56 and implemented in the MOLPRO pro¬ 
gram suite. 57 

In the first subsection of Section II, we pro¬ 
vide a quick overview of the range separated hy¬ 
brid + random phase approximation (RSH+RPA) 
method. Our derivation of the gradient of the 
RSH+RPA total energy is applicable for numerous 
variants of the long-range RPA correlation energy: 
each of these variants will be shortly discussed. 
The second subsection explains the construction 
of the Lagrangian for a generic RSH+RPA method 
and, in the third subsection, we present explicit 
working expressions for the gradients. Section III 
provides some practical details about our imple¬ 
mentation, followed by illustrative applications on 
geometry optimizations of small molecules and of 
charge transfer complexes. The paper is closed by 
conclusions and a short outlook of future develop¬ 
ments. We use atomic units throughout the whole 
paper. 

Theory 

RSH energy 

In the range-separated hybrid DFT framework one 
starts with a self-consistent independent-particle 
calculation by minimizing the energy 

£ rsh = min{(0|f + V ne + K!;|d>) +£g xc [n*]}, 

( 1 ) 

with respect to the molecular orbitals <)>,■( r) of the 
single determinant <f>. Here n ( p is the density asso¬ 
ciated with the single-determinant wave function, 
T is the kinetic energy operator, Vj ie is the nuclear 
attraction operator, V} r e is a long-range electron- 
electron repulsion operator constructed with the 
error function v l * e (r) = erf(/lr)/r, and £g xc [n<p] 
is the associated short-range Hartree-exchange- 
correlation density functional, written as: 

£hxcK] = JdrF(S(r)), (2) 
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where E, (r) is an array of quantities such as the 
density 


?7(r) - £ 4 0) <f)*(r ) fc(r) 
i 

= Y.D^x;(r)x,(r), (3) 

and other density-related parameters (spin density, 
reduced density gradients, etc ...) entering in the 
definition of the functional. <^4 will designate a 
component of ( resp. D(°)) is the RSH den¬ 

sity matrix in the molecular orbital basis (resp. in 
the atomic orbital basis). Throughout the paper, 
the indices p, q, r, s are used for general molecular 
orbitals, i, j, k for occupied molecular orbitals and 
a, b, c for virtual molecular orbitals; the indices p, 
v, p, (7, are used for atomic orbitals. 

It will be convenient to express the RSH energy 
with a fockian f = h + g sr + g lr , where the short- 
and long-range two-electron contributions are: 


8 
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sr 
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lr 

pq 



dUr) 

dd^pq 

lr - h{pr\qsf) , 


(4a) 


(4b) 


and the two-electron integrals follow the chemist’s 
notation. 

Using the above-defined quantities, the RSH en¬ 
ergy of Eq. (Eq. (1)) reads as: 

£ RS H = tr{d"»f}-ltr{dl°'g lr } 

+ (££>„]-trjd'V}), (5) 


where the last terms are the double-count correc¬ 
tion of the long-range Hartree-Fock energy 

— |tr jd(°)g lr j) on the one hand, and the sum of 
the short-range Hartree and exchange-correlation 
energies (age = £ , (j xc [no] — tr |d l0i g sr |), on the 
other hand. 


RPA energy 

Since the minimizing RSH wave function is a 
single-determinant approximation to the exact 
wave function, the long-range part of the RSH 
energy contains only Hartree and exchange terms. 
In principle, the exact ground-state energy can 
be obtained from the RSH energy by adding the 
long-range correlation energy E]!\ 

E = E rsh + E ] J, (6) 

provided the exact short-range exchange- 
correlation functional were known and we could 
solve exactly the long-range correlation problem. 
In this work we are interested in the random phase 
approximation (RPA) to the long-range correlation 
energy E'J. As it has been discussed in previous 
works , 43,53 several alternative RPA variants exist, 
and most of them can be expressed either in an 
adiabatic-connection formalism , 53 or can be re¬ 
formulated as ring approximations in the coupled 
cluster doubles (CCD) framework . 12 Among the 
numerous variants proposed in the CCD frame¬ 
work , 43 we focus our attention to the following 
ones: 

1. Direct RPA (dRPA-I) 

AdRPA-l = 5tr{‘K lrl TS RPA }, (7) 

with the dRPA amplitudes 1 T|( rpa satisfying 
the dRPA Riccati equations: 

RdRPA = (I + 1 T|| R p A ) 1 K h (I + 1 T|j R p A ) 

+ 1 TdRPA e + e 1 T|[ rpa = 0 (8) 

2. Direct RPA with SOSEX (SOSEX) 

^tvSOSEX — 2 tr { 1 ® 1 * * ^dRPA } • (9) 

where the amplitudes satisfy the previous 
dRPA Riccati equations, Eq. (Eq. ( 8 )). 

3. Approximate exchange RPA (RPAX2) 

£< Ii .rpax2 = 2 tr {* K 1 ' 1T RPAX j • (10) 

with the RPAX *T RPAX amplitudes satisfy- 
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ing the following Riccati-like equation: 

Rrpax = (1+ 'Trpax) 1 B lr (I+ lr lRp A x) 
+ lr lRPAX e + e * T RPAX = (11) 

This last variant has been suggested by Hes- 
selmann in a rCCD formalism. 28 An analo¬ 
gous approach can be derived in an adiabatic 
connection framework as well. 58 

The above equations are expressed in terms 
of symmetry adapted super-matrices for a closed 
shell system. The matrix elements of £ are £i a jb = 

f a hS,j — fij8 a j ), 1 K 1 ' and 1 B* 1 are combinations of 
two-electron integrals over spatial molecular or¬ 
bitals: 

l Kj' aJb = 2(ia\jb) h = 2K? aJb (12) 

= 2{ia\jbf-(ib\jaf = 2K* Jb -K% Jb . 

(13) 

RSH+RPA Lagrangian 

Since the RPA correlation energy is non- 
variational, we use the Lagrangian formalism to 
evaluate the RSH+RPA gradients. The Lagrange 
functional is constructed from the total RSH+RPA 
energy expression of Eq. (Eq. (6)) along with a set 
of undetermined Lagrange multipliers associated 
with the constraints that the parameters entering 
the energy expression must fulfill. In our case 
we have to consider three constraints which en¬ 
sure that (1) the orbitals are solutions of the RSH 
equations, (2) they remain orthogonal and (3) the 
amplitudes defining the RPA correlation energy 
are always solutions of the Riccati-like equations. 
The Lagrangian then reads: 

2^(C, T,z,x, A) = Ersh(C) + 4 r PA (C, T) 

+ tr {zf} + tr |x (C T SC — l) } 

+ tr{AR(C,T)}, (14) 

where C is the matrix of the MO coefficients ob¬ 
tained from the self-consistent RSH equations in 
the LCAO framework, T is the super-matrix of the 
RPA amplitudes in MOs, z is the set of multipli¬ 
ers associated with the Brillouin conditions, x is 
the multiplier matrix for the orthogonality condi¬ 


tions and A is the super-matrix related to the Ric- 
cati conditions R(C, T) . For the sake of notational 
simplicity, the the following the superscript "lr" 
will be omitted for quantities, where there is no 
risk of confusion (like RPA amplitude matrices). 
Explicit "lr/sr" labels will be kept only for cases 
where their use seemed indispensable for under¬ 
standing. 

The Lagrange multipliers A are obtained from 
the stationary condition of the Lagrangian with re¬ 
spect to the amplitudes T, which can be written 
for all the versions of RPA introduced in Section 
Eq. (5) as: 

Q(T)A + AQ(T) t = -P, (15) 

where the expression of the super-matrix Q(T) de¬ 
pends on the Riccati equation considered and the 
right-hand side P is a combination of two-electron 
integrals, having a form that depends on the corre¬ 
lation energy expression used for £|(|, A . 

After some lengthy algebra outlined in the Ap¬ 
pendix Table 4, one can show that the stationary 
conditions of the Lagrangian with respect to the 
orbital coefficients C (more precisely: with respect 
to the first-order variation of these coefficients, V) 
boil down to the following two equations: 

^© — © T + fz — zf 

< +4g lr [z]+4w sr [z]) = 0 (16a) 

/ ai 

k © + © T -|-@(z) -f-®(z) T = — 4x (16b) 

The first equation is a coupled perturbed RPA 
equation, which should be solved for z, the differ¬ 
ence density matrix in the virtual-occupied block; 
the second equation gives x, the energy-weighted 
density matrix, once z is known. The terms that 
depend neither on z nor on x are collected in the 
matrix ©, while 0(z) regroups terms that depend 
on z. 

As it is shown in the Appendix, there is a close 
correspondence between the procedures used to 
derive the long- and short-range contributions to 
Eq. (Eq. (16a)) and Eq. (Eq. (16b)) due to the ap¬ 
pearance of two analogous terms, denoted by g lr [z] 
and w sr [z] respectively. Furthermore, in the ex¬ 
pression of the Lagrange multiplier, x, defined by 
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Eq. (Eq. (16b)) one can identify contributions from 
the derivatives of both the RPA long-range corre¬ 
lation and the RSH reference energies, i.e. one can 
write for the occupied-occupied block: 

My = (*rpa) ( . —2(f)y. (17) 

where the second term is the usual expression of 
the x multiplier coming from the gradient of the 
reference energy. 

Analytical gradients 

Once the equations (Eq. (15)), (Eq. (16a)) and 
(Eq. (16b)) have been solved for the multipliers 
A, z and x, the Lagrangian is fully known. By 
the virtue of its variational property with respect 
to all of its parameters and since at its minimum 
it is equal to the RSH+RPA energy, the gradient is 
given by: 

£<*> =.Sf<'> = tr [ I) 1 Hj + tr | } + TU^p,' 

+ £ ( D2 + r2 ) w >''|P'7) lrW , 

[IVpG 

(18) 


usually in various types of gradient expressions, 
( D 2 ),v.ap = (l D( ° } + D(2) + Z ) Mv <a 

_I^i d (°) + d( 2 )+z)^dSS, 

( 20 ) 

and those contributions, which are specific to the 
RSH+RPA gradient expression: 

( r 2 )„v.c P = 

ia.jb 

+ 52 ClijCviCjaCjp (N ) ia j b . 

iajb 

( 21 ) 

The expressions of the super-matrices M and 
N depend on the Riccati equation and long-range 
correlation energy formula, specific to a given 
RPA variant. These quantities are coming from 
the factorization of the Lagrangian with respect to 
terms that depend on the orbital coefficients. The 
gradient expressions of the various RPA versions 
considered in this paper differ only in the details 
of the M and N super-matrices. For example, in 
the case of dRPA-I their expressions are: 


with (X) MV = LpqCup (\)p q Cj jV . and: 

( D V = E<w (<i ( 0 ) +<J ( 2 ) +z) M cfv 

= (d^+D^+z) . (19) 

As mentioned previously, the matrix of the 
Lagrange multipliers, x, can be identified as 
the energy-weighted one-particle density matrix, 
while z is the difference density matrix. The ma¬ 
trix d' 2 -* is defined in the Appendix Table 4. For 
the sake of clarity, in the four-index relaxed two- 
particle density appearing in equation (Eq. (18)) 
we have separated in the contributions that appear 


M 

N 


x 1 TdRPA + A + A 1 TdRPA + 1 TdRPA A 


+ 1 TdRPA A 'T d RPA (22) 

0. (23) 


The derivative of the density functional part with 
the corresponding double-counting term reads, on 
a real-space grid of points {A} which have the in¬ 
tegration weights (Ox : 
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*£?=E {e<»a +11,0)} 


X A 


dF 




, rr,v ^ /ed‘ 0) (.t) ed |2) (.t) ez(.v; 

+ LL®^^i + $A,A + $A,A 


+EE^(^+^ 


A AS 


:(•*) 

a- 


(24) 


Further details about the short-range DFT gra¬ 
dient contributions can be found in Appendix 
Eq. (33). Note that in the limiting case of van¬ 
ishing range-separation parameter, /A, one obtains 
the full-range DFA gradients, for /A —>■ °°, one gets 
the Hartree-Fock + RPA gradients, as given by 
Rekkedal et al., and by omitting completely the 
long-range RPA at any finite value of /A, one gets 
simply the analytical gradient expression of the 
RSH energy. 


Tests and applications 

The above equations have been implemented in 
the development version of MOLPRO. 57 Tak¬ 
ing advantage of some similarities in the struc¬ 
ture of the RSH+RPA and RSH+MP2 gradients, 
our implementation closely follows the flowchart 
of the MOLPRO MP2 gradient program 59 pre¬ 
viously adapted for the range-hybrid methods. 56 
Computational times are as expected, i.e. at the 
sr-LDA+lr-RPA level there is the same percent¬ 
age difference between the timings of the energy 
and gradient calculations as at the sr-LDA+lr-MP2 
level. There is also the same percentage additional 
cost when comparing an sr-LDA+lr-MP2 calcula¬ 
tion to an sr-LDA+lr-RPA calculation of either the 
energy or the gradients. The scaling of the gradi¬ 
ent calculation follows that of an energy calcula¬ 
tion (in the present implementation A 6 ). However 
the algorithm can take easily advantage of future 
density-fitting implementations, which can be as 
low as N 4 for the dRPA-I method and N 5 for other 
RPA variants. 

Lor all the data presented here, we calculate the 


mean absolute error (MAE), defined as ^ ]T ( \a t — 
^y f |, and the mean percentage absolute error de¬ 
fined as 

Table Table 1 and Ligure Ligure 1 show bond 
lengths and angles obtained from geometry opti¬ 
mizations of a set of small molecules (H 0 , HE, 
H 2 0, HOF, HNC, NH 3 , N 2 H 2 , HNO, C 2 H 2 , HCN, 
C 2 H 4 , CH 4 , N 2 , CH 2 0,"cO, C02, F 2 ) at the 
sr-LDA+lr-MP2, sr-LDA+lr-dRPA-I, sr-LDA+lr- 
SOSEX as well as sr-LDA+lr-RPAX2 levels. Op¬ 
timizations were performed using the program 
GADGET 60 interfaced with MOLPRO, with the 
aug-cc-pVQZ basis set, although, as expected in 
an RSH framework, the convergence of the results 
is fast with respect to the basis set. The calcu¬ 
lations were considered to have reached conver¬ 
gence when all gradients components were under 
0.0003 a.u. Our results are compared to the opti¬ 
mized geometries, recently published by Re kk edal 
et al, 51 at the MP2, dRPA-I and SOSEX lev¬ 
els without range-separation, and perfectly repro¬ 
duced by our present implementation. The refer¬ 
ence geometries have been taken from the work of 
Pawlowski et al. 61 and are obtained from experi¬ 
mental rotational constants and vibration-rotation 
interaction constants computed at the CCSD(T) 
level with the cc-pVQZ basis. The deviations from 
the reference bond lengths are usually less than 
4 pm, except the cases of F—F and O—F bond 
lengths, where the error can be as large as 7 pm. 
In the calculations using RHF orbitals, the mean 
absolute error (MAE) of the MP2 bond lengths is 
0.476 pm while dRPA-I and SOSEX have mean 
absolute errors of 1.589 pm and 2.077 pm, re¬ 
spectively; the range-hybrid calculations, on the 
other hand, all yield similar MAE (1.276 pm for 
sr-LDA+lr-MP2, 1.347 pm for sr-LDA+lr-dRPA-I, 
1.459 pm for sr-LDA+lr-SOSEX and 1.388 pm for 
sr-LDA+lr-RPAX2). Among all the methods pre¬ 
sented in Table Table 1, it is the full-range MP2 
which gives the best results, at least at this rel¬ 
atively high basis set level. A comparison with 
the MAE values reported by Burow at al. 52 indi¬ 
cates that simple PBE and PBE0 functionals have 
a performance which is not very far from that of 
full-range MP2. Their dRPA-I calculations using 
PBE orbitals is slightly better than MP2. It is un¬ 
deniable that the RSH+RPA version tested here, 
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based on short-range LDA, have a worse perfor¬ 
mance than these above-mentioned methods. The 
range-hybrid sr-LDA+lr-dRPA-I and sr-LDA+lr- 
SOSEX calculations yield globally better agree¬ 
ment with the reference values then analogous cal¬ 
culations using RHF orbitals. However, it is quite 
clear that mainly the X—H bond lengths are im¬ 
proved in the range-hybrid calculations, while the 
bond lengths between first-row atoms remain es¬ 
sentially of similar quality as in RHF+RPA calcu¬ 
lations. The global agreement is essentially the 
same between the four range-separated methods 
(MP2 and dRPA-I, SOSEX, RPAX2) indicating 
that in these simple cases the bond length cor¬ 
rections, which can be attributed to higher than 
second order Mpller-Plesset effects are negligible. 
This is supported by the fact that different ways of 
including exchange diagrams in the correlation en¬ 
ergy calculations, namely following either SOSEX 
or RPAX2, does not improve the results signifi¬ 
cantly. The angles are somewhat farther on from 
the reference in the range-hybrid cases than in the 
calculations with RHF orbitals, but the differences 
are only in the order of one degree. 



Figure 1: Deviation of bond distances (pm) of 
small molecules with respect to the reference. 61 

In order to test the gradients for a class of inter- 
molecular interactions, we present interaction en¬ 
ergies (Table Table 2 and Figure Figure 3) as well 
as inter-monomer distances (Table Table 3 and 
Figure Figure 4) resulting from the geometry op¬ 
timizations of binary systems in the CT7 (charge 
transfer) ensemble of intermolecular complexes 62 
at the sr-FDA+lr-dRPA-I, sr-FDA+lr-SOSEX and 
sr-FDA+lr-RPAX2 levels. For an illustration of 
the relative orientation of the monomers, see Fig¬ 
ure Figure 2. Our results are compared to those 


from Chabbal et al. 56 obtained by geometry op¬ 
timizations at the MP2 and sr-FDA+lr-MP2 levels 
and to reference values given by the group of Truh- 
lar. 62,63 The optimizations were conducted with 
GADGET program using the aug-cc-pVTZ basis 
set, without counterpoise correction. It is believed 
that for range-hybrid calculations the effect of the 
basis set superposition error is small, even at this 
relatively modest basis set level. 

These charge transfer complexes are usually 
problematic for plain DFT methods due to the siz¬ 
able delocalization error of the common function¬ 
als, as it has been well-known for a long time. 64 
Therefore it is expected that the RSH determi¬ 
nant is going to be a reasonable reference state 
for the long-range RPA correlation calculations to 
take into account the Fondon dispersion interac¬ 
tions stabilizing these complexes. The study from 
Chabbal et al. has already shown that the range- 
separated MP2 approach improves the results with 
respect to full-range MP2 for charge transfer com¬ 
plexes. Our results demonstrate that the three lr- 
RPA variants tested yield a general improvement 
over RSH+MP2 calculations, with mean absolute 
errors around 0.30 kcal.mol 1 for the interaction 
energies and around 2.5 pm for the inter-monomer 
distances. 

The percentage errors of all methods for both the 
interaction energy and the inter-monomer distance 
is rather high in the case of the NH 3 ■ • • F 2 dimer 
(c.f. for example the visual abstract which shows 
the percentage deviation of the inter-monomer dis¬ 
tances). This observation can be attributed to the 
relatively small magnitude of the reference val¬ 
ues, especially in the case of inter-monomer dis¬ 
tances. While the H 9 0 ■ C1F interaction energy 
is underestimated by sr-FDA+lr-dRPA-I and well 
recovered by sr-FDA+lr-SOSEX and sr-FDA+lr- 
RPAX2, the inter-monomer distance after geom¬ 
etry optimization is consistently less good using 
any of the range-hybrid RPA method, with er¬ 
rors around 4.5 pm. On the contrary, while the 
inter-monomer distance of the NH 3 • ■ • Cl 9 dimer 
is close to the reference for all the range-hybrid 
RPA methods, the interaction energies show the 
largest deviations (around 0.6 kcal.mol -1 ) in the 
whole set of systems. We are going to attempt 
a rationalization of these observations in the next 
paragraph. 
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Table 1: Deviation with respect to the reference of the bond lengths r e (pm) and angles a (degrees) of 
small molecules obtained after geometry optimizations, "from Ref. 51 




sr-LDA+ 

lr-MP2 

sr-LDA+ 

lr dRPA I 

sr-LDA+ 

lr-SOSEX 

sr-LDA+ 

lr-RPAX2 

MP2“ 

dRPA-I a 

SOSEX" 

Ref. 61 

H 2 

H H 

0.572 

0.783 

0.581 

0.616 

-0.543 

-0.814 

-0.794 

74.149 

HF 

F-H 

0.323 

0.394 

0.427 

0.296 

-0.077 

-1.067 

-1.491 

91.688 

h 2 o 

O-H 

-0.120 

-0.099 

-0.202 

-0.192 

-0.182 

-1.139 

-1.504 

95.790 

HOF 

O-H 

-0.203 

-0.091 

-0.454 

-0.269 

-0.399 

-1.510 

-1.908 

96.862 

HNC 

N-H 

0.219 

0.325 

0.260 

0.150 

-0.099 

-1.014 

-1.143 

99.489 

nh 3 

N-H 

-0.394 

-0.383 

-0.480 

-0.496 

-0.386 

-1.074 

-1.334 

101.139 

n 2 h 2 

N-H 

-0.185 

0.144 

0.470 

-0.279 

-0.403 

-1.285 

-1.555 

102.883 

HNO 

N-H 

-0.054 

-0.233 

-0.365 

-0.172 

-0.557 

-1.435 

-1.765 

105.199 

c 2 h 2 

C-H 

0.156 

0.185 

0.094 

0.092 

-0.253 

-0.811 

-0.865 

106.166 

HCN 

C-H 

0.128 

0.175 

0.058 

0.055 

-0.331 

-0.925 

-0.940 

106.528 

c 2 h 4 

C-H 

0.016 

0.332 

0.288 

-0.056 

-0.372 

-0.920 

-0.942 

108.068 

ch 4 

C-H 

-0.131 

-0.096 

-0.241 

-0.227 

-0.431 

-0.840 

-0.873 

108.588 

n 2 

N-N 

-1.644 

-1.751 

-1.865 

-1.847 

1.010 

-1.552 

-2.209 

109.773 

ch 2 o 

C-H 

0.148 

0.151 

0.071 

0.105 

-0.404 

-1.008 

-0.977 

110.072 

CO 

C-O 

-1.392 

-1.454 

-1.574 

-1.556 

0.370 

-1.392 

-1.914 

112.836 

HCN 

C-N 

-1.781 

-1.921 

-2.095 

-2.062 

0.711 

-1.593 

-2.178 

115.336 

co 2 

C-O 

-1.161 

-1.235 

-1.335 

-1.322 

0.353 

-1.364 

-1.872 

116.006 

HNC 

C-N 

-1.583 

-1.698 

-1.782 

-1.796 

0.160 

-1.413 

-1.882 

116.875 

c 2 h 2 

c-c 

-1.649 

-1.830 

-2.002 

-1.964 

0.243 

-1.464 

-1.900 

120.356 

ch 2 o 

C-O 

-1.641 

-1.908 

-1.925 

-1.797 

0.105 

-1.632 

-2.224 

120.465 

HNO 

O-N 

-2.778 

-2.811 

-2.868 

-2.900 

0.611 

-2.436 

-3.355 

120.859 

n 2 h 2 

N-N 

-2.780 

-3.019 

-3.095 

-2.973 

0.291 

-2.188 

-3.008 

124.575 

c 2 h 4 

C-C 

-1.901 

-2.493 

-2.475 

-2.077 

-0.464 

-1.478 

-1.811 

133.074 

f 2 

F-F 

-5.730 

-5.534 

-5.978 

-5.941 

-1.737 

-4.998 

-7.208 

141.268 

HOF 

O-F 

-5.214 

-4.634 

-5.479 

-5.460 

-1.412 

-4.376 

-6.270 

143.447 

HOF 

H-O-F 

2.214 

-1.807 

2.407 

2.389 

0.138 

1.453 

2.178 

97.860 

h 2 o 

H O-H 

1.886 

2.089 

2.199 

2.186 

-0.249 

0.472 

1.065 

104.400 

n 2 h 2 

H-N-N 

1.605 

1.892 

1.296 

1.666 

-0.435 

0.587 

1.075 

106.340 

nh 3 

H-N-H 

1.757 

1.900 

1.982 

2.032 

-0.527 

-0.241 

0.320 

107.170 

HNO 

H-N-O 

1.008 

1.151 

1.217 

0.946 

-0.495 

0.296 

0.671 

108.260 

c 2 h 4 

C-C H 

0.164 

-0.016 

0.707 

0.176 

-0.063 

0.128 

0.127 

121.400 

ch 2 o 

H-C-O 

0.098 

0.401 

0.180 

0.115 

0.146 

0.264 

0.283 

121.630 

r e MAE 

1.276 

1.347 

1.459 

1.388 

0.476 

1.589 

2.077 


r e M%AE 

1.034 

1.103 

1.191 

1.128 

0.415 

1.373 

1.777 


a MAE 


1.247 

1.322 

1.427 

1.359 

0.293 

0.492 

0.817 


a M%AE 

1.195 

1.258 

1.356 

1.302 

0.273 

0.473 

0.787 
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Table 2: Interaction energies (kcal.mol 1 ) of the CT7 dimers after geometry optimization without coun¬ 
terpoise correction. "results from. 56 


MP2 a 

sr-LDA+ lr-MP2 fl 

sr-LDA+ 

sr-LDA+ 

sr-LDA+ 

Ref. 62 



Ir-dRPA-I 

lr-SOSEX 

lr-RPAX2 



c 2 h 4 - 

•f 2 

1.56 

1.16 

0.95 

0.80 

0.82 

1.06 

nh 3 - 

■ -f 2 

1.99 

1.71 

1.31 

1.32 

1.33 

1.81 

c 2 h 2 -- 

•C1F 

4.89 

4.36 

3.41 

3.53 

3.55 

3.81 

HCN- • 

•C1F 

5.72 

5.81 

4.96 

5.06 

5.08 

4.86 

nh 3 - 

•ci 2 

5.59 

5.19 

4.25 

4.26 

4.29 

4.88 

h 2 o-- 

•OF 

6.00 

6.29 

5.01 

5.41 

5.43 

5.36 

nh 3 -- 

•OF 

11.95 

12.10 

10.62 

10.66 

10.72 

10.62 

MAE 


0.76 

0.63 

0.30 

0.28 

0.28 


M%AE 


20.31 

12.37 

10.05 

10.98 

10.70 



Table 3: Inter-monomers distances (pm) of CT7 dimers after geometry optimization without counterpoise 
correction, "results from. 56 


MP2 a 

sr-LDA+ lr-MP2 fl 

sr-LDA+ 

sr-LDA+ 

sr-LDA+ 

Ref 62.63 



Ir-dRPA-I 

lr-SOSEX 

lr-RPAX2 



c 2 h 4 - 

•f 2 

291.2 

301.4 

305.0 

305.7 

305.9 

305.3 

nh 3 - 

• f 2 

264.6 

267.8 

273.3 

275.2 

274.9 

269.6 

c 2 h 2 -- 

•C1F 

280.0 

282.7 

289.1 

289.6 

288.6 

287.6 

HCN- • 

•C1F 

254.8 

253.2 

263.1 

258.7 

258.7 

260.9 

nh 3 -- 

•Cl 2 

260.3 

264.4 

271.8 

270.2 

270.0 

268.8 

h 2 o- • 

•C1F 

251.2 

247.6 

251.0 

251.5 

251.8 

255.7 

nh 3 -- 

•C1F 

223.4 

224.1 

231.5 

228.4 

228.3 

230.2 

MAE 


7.5 

5.3 

2.4 

2.5 

2.3 


M%AE 


2.8 

2.0 

0.9 

1.0 

0.9 
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Q 


> 


0=0 


NH 3 -F 2 c 2 h 4 -f 2 


> 


nh 3 --ci 2 c 2 h 2 -cif 




HCN--C1F NHv-CIF 


H 2 0- -C1F 


Figure 2: Charge transfer complexes of the CT7 
set. 62 



Figure 5: Effect of the dimerization on the bond 
lengths (pm) of the monomers of the CT7 dimers. 
The reference is from. 62,63 



C 2 H4-F 2 NH 3 -F 2 C 2 H 2 -C1F HCN-C1F NH 3 -C1 2 H 2 0-C1F NH 3 -C1F 


Figure 3: Deviation of the interaction energies 
(kcal.mol -1 ) of the CT7 dimers after geometry 
optimization without counterpoise correction with 
respect to the reference. 62 



Figure 4: Deviation of the inter-monomers dis¬ 
tances (pm) of CT7 dimers after geometry opti¬ 
mization without counterpoise correction with re¬ 
spect to the reference. 62,63 


We propose on Table Table 4 and Figure Fig¬ 
ure 5 an analysis of the bond lengths obtained after 
geometry optimizations via the difference Adi m = 
v dimer—Anono between the bond lengths optimized 
in the dimer (rdi mer ) and the independently opti¬ 
mized bond lengths in the monomers (r mono ). This 
quantity measures the effect of the dimerization 
on the geometry of the monomers. The rdi mer and 
Adim values resulting from sr-LDA+lr-dRPA-I, sr- 
LDA+lr-SOSEX and sr-LDA+lr-RPAX2 calcula¬ 
tions and for the reference geometries from Zhao 
and Truhlar 62,63 obtained at the MC-QCISD/3 
level, are shown in Table Table 4. The bonds 
are collected in different groups: C—X, X—H, 
X—X (heteroatomic and homoatomic). We ob¬ 
serve that the bonds are generally slightly less de¬ 
formed by the intermolecular interactions in the 
case of the sr-LDA+lr-dRPA-I calculations than 
in reference geometries, and are better recovered 
at the sr-LDA+lr-SOSEX and sr-LDA+lr-RPAX2 
levels (see Figure Figure 5). Both bonds involved 
in the dimer NH 3 • • • F 2 reflect the standard behav¬ 
ior of their groups. This confirms that this com¬ 
plex is not an inherently problematic case in range- 
hybrid RPA. The NH 3 ■ ■ • Cl 9 complex shows in¬ 
ternal bond lengths far off as compared to the stan¬ 
dard behavior. The bonds are more deformed in 
the RSH-RPA calculations in comparison to the 
reference geometry, which explains a large un¬ 
derestimation of the interaction energy in spite 
of the good inter-monomer distances. We see 
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that the bonds involved in the FP,0 ■ • • C1F dimer, 
which showed bad inter-monomer distances for all 
RSH-RPAs, are described much better by the sr- 
LDA+lr-SOSEX and sr-LDA+lr-RPAX2 methods 
than by the sr-LDA+lr-dRPA-I: this could explain 
the improvement previously mentioned for the in¬ 
teraction energy of this dimer. 
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Table 4: Data concerning the bond lengths (pm) of the monomers of the CT7 dimers, /'dimer are the bond lengths optimized in the dimer. Adj m is the 
effect of the dimerization on the bond lengths (Adi m > 0 shows a bond that is longer in the dimer), AA = A^ hod — Ajj[ T1 is the compared effect of the 
dimerization between a given method’s level and the MC-QCISD/3 level. 



sr-LDA+lr-dRPA-I 

F dimer ^dim AA 

sr-LDA+lr-SOSEX 

r dimer Adim AA 

sr-LDA+lr-RPAX2 

f dimer Adi m AA 

Ref. 62 ’ 63 

f dimer Adi m 

HCN- • C1F 

C-N 

113.363 

-0.173 

0.016 

113.171 

-0.166 

0.023 

113.206 

-0.166 

0.023 

115.621 

-0.189 

C 2 H 2 ---C1F 

u 

1 

u 

118.788 

0.179 

-0.039 

118.613 

0.182 

-0.036 

118.656 

0.181 

-0.037 

121.058 

0.218 

c 2 h 4 --f 2 

u 

1 

u 

131.237 

0.368 

0.306 

131.025 

0.042 

-0.020 

131.050 

0.029 

-0.033 

133.678 

0.062 

c 2 h 2 -cif 

C-H 

106.321 

-0.051 

-0.123 

106.362 

0.090 

0.017 

106.400 

0.120 

0.047 

106.668 

0.073 

nh 3 -f 2 

N-H 

100.465 

-0.375 

-0.405 

100.741 

0.025 

-0.005 

100.753 

0.025 

-0.005 

101.524 

0.030 

c 2 h 4 ---f 2 

C-H 

108.111 

-0.585 

-0.576 

108.045 

-0.006 

0.003 

108.054 

-0.008 

0.000 

108.437 

-0.009 

h 2 o-cif 

O-H 

95.356 

-0.452 

-0.570 

95.853 

0.149 

0.031 

95.861 

0.147 

0.028 

96.301 

0.118 

HCN •••C1F 

C-H 

106.658 

-0.107 

-0.135 

106.692 

0.094 

0.066 

106.699 

0.094 

0.066 

106.930 

0.028 

NHy-ClF 

N-H 

101.121 

0.280 

0.339 

100.870 

0.154 

0.212 

100.882 

0.153 

0.212 

101.435 

-0.058 

nh 3 ---ci 2 

N-H 

100.931 

0.090 

0.058 

100.797 

0.081 

0.050 

100.810 

0.081 

0.050 

101.525 

0.031 

HCN-C1F 

Cl—F 

162.386 

0.992 

-0.384 

162.263 

1.608 

0.232 

162.320 

1.612 

0.236 

165.635 

1.376 

NHyCIF 

Cl—F 

166.718 

5.325 

-0.511 

166.796 

6.140 

0.304 

166.901 

6.193 

0.357 

170.095 

5.836 

C 2 H 2 -C1F 

Cl—F 

162.159 

0.766 

-0.652 

162.044 

1.389 

-0.029 

162.131 

1.422 

0.004 

165.677 

1.418 

h 2 o-cif 

Cl—F 

162.956 

1.563 

0.121 

162.372 

1.717 

0.275 

162.419 

1.711 

0.269 

165.701 

1.442 

NH 3 • Fj 

F-F 

136.602 

0.740 

-0.384 

136.226 

0.821 

-0.303 

136.295 

0.833 

-0.291 

142.547 

1.124 

c 2 h 4 ---f 2 

F-F 

136.020 

0.159 

-0.414 

135.849 

0.443 

-0.130 

135.900 

0.438 

-0.136 

141.996 

0.573 

nh 3 ---ci 2 

Cl-Cl 

199.898 

2.303 

-0.193 

199.388 

2.363 

-0.133 

199.526 

2.395 

-0.101 

203.628 

2.496 











Conclusions and outlook 

The RSH+RPA analytical energy gradients have 
been derived using the Lagrangian formulation 
and implemented in the development version of 
the MOLPRO quantum chemical program pack¬ 
age. Although the working expressions have been 
obtained for all of the main categories of the RPA 
correlation energy, the present work reports the nu¬ 
merical implementation only for the direct RPA, 
the SOSEX and the RPAX2 variants. These results 
show a significant improvement with respect to 
range-separated MP2 calculations in the descrip¬ 
tion of both the energetics and the structure of 
charge transfer complexes, where intermolecular 
interactions play an important role even in the ge¬ 
ometry of the constituents of the complexes. 

The present numerical implementation provides 
mainly reference data for relatively small systems. 
Admittedly, the computational efficiency of an 
orbital-based algorithm, used here, is quite lim¬ 
ited. However, generalizations for density-fitting, 
resolution-of-identity and even Cholesky decom¬ 
position algorithms (outlined e.g.in Ref. 12 ) seem 
to be rather straightforward and will be the sub¬ 
ject of future work. Another extension of the 
present work consists in the computational realiza¬ 
tion of the exchange-including rCCD/RPAx cor¬ 
relation energy expressions, namely the SOI and 
S02 variants, 43 which have shown the best quali¬ 
tative performance in range-hybrid calculations of 
intermolecular interaction energies. 

It is quite clear from our past experience, that 
for simple intermolecular interact action energies 
the short-range functional has a relatively minor 
influence on the quality of the results. The sit¬ 
uation seems to be different as far as we would 
like to reproduce bond lengths and angles, and 
the use of short-range GGA (e.g. sr -PBE) function¬ 
als is mandatory to improve these results. Note 
that in geometry optimizations, in addition to 
the sr exchange-correlation functionals and poten¬ 
tials one needs sr-PBE kernels (second functional 
derivatives of the sr-PBE functionals). Work in 
this direction is in progress and we hope to test 
this hypothesis in the near future. 

As a by-product of the analytical force imple¬ 
mentation, the non-relaxed and relaxed density 
matrices at the RSH-RPA levels are available. 


They will be exploited for the analysis of the cor¬ 
relation effects on one-electron properties, like 
charge densities and their multipole moments. A 
study in this direction is in preparation. 
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Stationary conditions with re¬ 
spect to orbital coefficients 

Let us parameterize the variation of orbital coeffi¬ 
cients at first order by a unitary rotation matrix V 
as C <— C + CV. The stationary conditions for the 
Lagrangian can be written as 


1 dS£ 

2 ~dV 


= 0 


v=o 


(25) 


The factor of \ is inserted to compensate a fac¬ 
tor of 2 appearing in the upcoming derivations for 
reasons of symmetry. 

It will prove to be convenient to rewrite the La¬ 
grangian of Eq. (Eq. (14)) by factorizing the terms 
which depend on the orbital coefficients. This 
can be achieved by separating in Er PA (C. T) + 
tr{AR(C,T)} the terms depending on the fock- 
ian from those depending on the two-electron in¬ 
tegrals, leading to: 


j£?(C,T,z,x,A) = tr{df}+B DC 

+ tr{MK}+tr{NK'} 

+ tr {x (C T SC — l) } , (26) 

where the super-matrices M and N gather all the 
elements that are multiplied by the integrals K and 
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K 7 , respectively. The particular forms of these 
super-matrices thus depend on the Riccati equa¬ 
tion and on the energy expression corresponding 
to the given RPA variant chosen for the long-range 
correlation energy. We defined the relaxed den¬ 
sity d = +d (2 ^ Tz with the matrix d* 2 ), whose 

blocks are: 

(d (2) ) y = -{T, - {A,T} y 

(d <2) ) o6 = {T.A} oi ,+ {A,T}„ i (27) 

(d< 2 >) — 0 

V / ai 

In the above equations we use a specific notation 
for the "failed traces", that is to say for the partial 
summations leading to a result which still depends 
on two of the four indexes that compose the super¬ 
indexes of the super-matrices: 


{X- Y}// — V, Xjg k c Ykc,ia 

(28) 

kc,a 


{X,Y} fl z, = V, X ia .kcYkc. t b 

(29) 

kc,i 



The derivatives of all the terms in Eq. (Eq. (26)) 
with respect to a change in the orbital coefficient 
are fairly lengthy, therefore only some of the ele¬ 
ments are given here. From the derivation of the 
trace of the two-electron integrals with the super¬ 
matrix M will emerge contractions of the form 
{k,m} (7 and {K,M} ai , as well as generalizations 
of the form {K, M). and {K, M) . where the 
super-matrix K is constructed from K but does not 
respect it’s iajb structure (the same quantities are 
derived for the trace with N). All those terms are 
grouped in the matrix ©. 

Two-electron fockian and double¬ 
count derivatives 

The derivatives of the long- and short-range 
parts of the "fockian plus double-count" terms 
show some interesting analogies. The deriva¬ 
tive of the long-range two-electron contribution 

tr|dg lr d (0) | yields a g lr d (0) d, and, by an 


interchange property of the indexes involved in 
the summations, g lr |d| d i0i . The derivation of 
the double-count correction will cancel out the 
glr d ( °) d(°) in the "interchanged" term, so that 
we finally obtain: 



In a similar, but less obvious manner, the deriva¬ 
tive of the short-range contribution tr{dg sr } gives 
g sr d and, by a comparable phenomenon, a new ob¬ 
ject that we call w sr [d] d (0) . Using the relationship: 


1 d^Hxc N „sr j(0) 

2 3\ v=0 8 ’ 


(31) 


(see Appendix Eq. (33)) we see that the double¬ 
count term behaves as follow: 


= g sr d(°) - g sr d(°) - W sr d (0) d (0) , 

(32) 

and cancels out the w sr d (0) d (0) from the 
"interchanged" term, much like previously the 
long-range double-count term cancelled out 
glr d (0 i d :0) from the long-range "interchanged" 
term. This leads to the result: 


z ° v v=o 

= g sr d + w sr d (2) +z d (0) (33) 

The relationships (Eq. (30)) and (Eq. (33)) are 
similar to each other and their terms enter in the 
definition of © and 0(z). 
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Elaboration on the sr-DFT terms 


In order to derive some of the quantities needed for 
the short-range treatment of the gradient, we write 
the Hartree-exchange-correlation functional as (cf. 
Eq. (Eq. (2))): 

£hxcM = fdrF(Z(r)), (34) 

where £ is an array of quantities that enter in 
the definition of the functional, i.e. = 

{n a ,V/ 7 a ,V/ 7 ^,...}. For all such objects, that do 
not include explicitly virtual orbitals, we can write 
in a most general way: 

& = [I ^n A («0) = FI E ( d pq QpQq) » 

n n pq 

(35) 

where ■%, are functions that are different for ev¬ 
ery quantity £ 4 . 

Initially, with these notations at hand, we are 
going to show that the derivative of the func¬ 
tional with respect to the orbital rotation param¬ 
eters is related to the two-electron part of the 
short-range fockian, i.e. we are to going to prove 
Eq. (Eq. (31)). The derivative of the functional 
with respect to the orbital rotation parameters 
reads: 


I^Hxc 

2 dV rs 

(36) 


1 y f r dF dU 

, =0 2 L*} 1 dU dv rs 


and the two-electron part of the short-range fock¬ 
ian is: 


8 


sr 

rs 


L 


3F d£, 

a4? 


(37) 


The derivative of a quantity £4 with respect to 
a change of the orbital coefficients, appearing in 
Eq. (Eq. (36)), is 



while the derivative with respect to d^ 0 ), which oc¬ 
cur in Eq. (Eq. (37)) is 


dU_ 

d4f 


i n^i 


(39) 


The above results allow us to write: 



Subsequently, we will show that the derivative 
of the trace of the short-range fockian, found in 
Eq. (Eq. (33)), is indeed the sum of terms g sr d and 
w sr [d] d :0: '. The derivative reads: 


r 1 f , 1 

h*J d ^ A dV rs Vrs =Qp q d d { pf Pq 


V „=0 

(41) 

Inspection of the element -^y d pq occurring 

ddpq 

in both terms of Eq. (Eq. (41)) reveals that it is a 
sum of objects E,a where one occurence of d^ 0) has 
been eliminated by the derivation and replaced by 
the density matrix d. We call this element it’s 
expression is: 


+E 4 /*• 


3F 3 


3%a 3V r , 


I^lu 


.pq 


dd 


( 0 )' 


pq 


pq 
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Eq. (Eq. (18)) is derived as: 


4? = I dp, = L^(»d)n & (no), 

pq odpq i riy^i 

(42) 

The first term of Eq. (Eq. (41)) involves only el¬ 
ements of the form "U rp d$" (see Eq. (Eq. (38))) 

that enter in the composition of w sr [d] d$ ■ The 
derivation of the second term in Eq. (Eq. (41)) is 
more involved. Two types of terms will arise from 
the derivative with respect to V of the element 
, as follows: 


psr(.v) _ F sr(x) 

C DFT — ^Hxc 


+ tr 






dUr) 

dd { 3 


(d( 2 )+ 

(45) 


With a quadrature on a grid of the real space 
of points {A} of weights (Ox , this leads to 
Eq. Eq. (24) of the main text. 


dV rs 


v„=o 


I^(nd)n^N) 


dV r , 


n^i 


y d^(n d ) 

h dv, 


\\^ A [n 0 ) 

rs V„= o n & 


v rs =o 


+E^(«d)I 

i j 


d^f A (n 0 ) 

dV rs 


n^i 

Vrs=O n £j 

(43) 


The first type of terms involves the derivative 
with respect to a rotation of the orbital coefficients 
of a && that has been "contaminated" by d: they 
raise elements of the kind "U rp d ps ” that contribute 
to g sr d; the second type of terms involves the 
derivative with respect to V of an "original" 3£*= A 
and will see the emergence of elements "U rp d$" 
that compose w sr [d] d (0 '. 

After all the derivations described above have 
been carried out, the full definition of w sr [d] is: 


+m (^a) n & («o) 

(44) 


With the notations introduced here, the deriva¬ 
tive of the short-range fockian term with re¬ 
spect to the atomic coordinates that appear in 
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